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successfully detect a terrestrial planet angularly close to its parent star: photon 
noise from diffracted star light, and speckle noise from star light scattered by 



ABSTRACT 



High-contrast imaging from space must overcome two major noise sources to 



instrumentally-generated wavefront perturbation. Coronagraphs tackle only the 
photon noise contribution by reducing diffracted star light at the location of a 
planet. Speckle noise should be addressed with adaptative-optics systems. Fol- 
lowing the tracks of Malbet, Yu, & Shao (1995), we develop in this paper two 
analytical methods for wavefront sensing and control that aims at creating dark 
holes, i.e. areas of the image plane cleared out of speckles, assuming an ideal 
coronagraph and small aberrations. The first method, speckle field nulling, is 
a fast FFT-based algorithm that requires the deformable-mirror influence func- 
tions to have identical shapes. The second method, speckle energy minimization, 
is more general and provides the optimal deformable mirror shape via matrix 
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inversion. With a, N xN deformable mirror, the size of matrix to be inverted 
is either N 2 x N 2 in the general case, or only iV x N if influence functions can 
be written as the tensor product of two one-dimensional functions. Moreover, 
speckle energy minimization makes it possible to trade off some of the dark hole 
area against an improved contrast. For both methods, complex wavefront aber- 
rations (amplitude and phase) are measured using just three images taken with 
the science camera (no dedicated wavefront sensing channel is used), therefore 
there are no non-common path errors. We assess the theoretical performance of 
both methods with numerical simulations including realistic speckle noise and 
experimental influence functions. We find that these speckle nulling techniques 
should be able to improve the contrast by several orders of magnitude. 

Subject headings: Instrumentation: adaptive optics — techniques: high angular 
resolution — planetary systems 

1. Introduction 

The field of extrasolar planet research has recently made a leap forward with the direct 
detection of extrasolar giant planets (EGPs). Using Spitzer Space Telescope, Charbonneau 
et al. (2005) and Deming et al. (2005) have detected infrared photons from two transit- 
ing planets, TrES-lb and HD209458b, respectively. Chauvin et al. (2004, 2005) have re- 
ported the infrared imaging of an EGP orbiting the nearby young brown dwarf 2M1207 with 
VLT/NACO, whereas Neuhauser et al. (2005) have collected evidence for an EGP companion 
to the T-Tauri star GQ Lup using VLT/NACO as well. 

Although there are claims that the direct detection of terrestrial planets could be per- 
formed from the ground with - yet to come - extremely large telescopes (Angel 2003; Chelli 
2005), it is widely believed that success will be more likely in space. Direct detection is the 
key to spectroscopy of planetary atmospheres and discovery of biomarkers, namely indirect 
evidence of life developed at the planetary scale (e.g. Des Marais et al. 2002). 

Both NASA and ESA have space mission studies well underway to achieve this task. 
Darwin, the European mission to be launched in 2015, will be a thermal infrared nulling 
interferometer with three 3.5-m free-flying telescopes (Karlsson et al. 2004). Terrestrial 
Planet Finder, the American counterpart, will feature two missions: a 8x3.5 m-monolithic 
visible telescope equipped with a coronagraph (TPF-C) to be launched in 2015, and an 
analog to Darwin (TPF-I) to be launched in the 2015-2019 range (Coulter 2004). 

The direct detection of the photons emitted by a terrestrial planet is made very chal- 
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lenging by the angular proximity of the parent star, and by the very high contrast (i.e. 
luminosity ratio) between the planet and its star: about 10 6 in the thermal infrared and 
about 10 10 in the visible. Both wavelength ranges have their scientific merits and technical 
difficulties, and both of them are thought to be necessary for an unambiguous detection of 
habitability and signs of life (e.g. Des Marais et al. 2002). In this paper, we deal with the 
visible range only. 

In the visible, planet detection faces two fundamental noise sources: (i) quantum noise of 
the diffracted star light, and (ii) speckle noise due to the scattering of the star light by optical 
defects. Labeyrie (1995) proposed a technique based on dark speckles to overcome speckle 
noise: random fluctuations of the atmosphere cause the speckles to interfere destructively 
and disappear at certain locations in the image, thus creating localized dark spots suitable 
for planet detection. The statistical analysis of a large number of images then reveals the 
planet as a spot persistently brighter than the background. 

Malbet, Yu, & Shao (1995) proposed to use a deformable mirror (DM) instead of the 
atmosphere to make speckles interfere destructively in a targeted region of the image called 
search area or dark hole (DH or H). Following the tracks of these authors, this paper dis- 
cusses methods to reduce the speckle noise below the planet level by using a DM and an 
ideal coronagraph. However, unlike Malbet, Yu, & Shao (1995), we propose non-iterative 
algorithms, in order to limit the number of long exposures needed for terrestrial planet 
detection. We will refer to these methods as speckle nulling techniques, as Trauger et al. 
(2004) call them. Technical aspects of this work are inspired by the High Contrast Imag- 
ing Testbed (HCIT; Trauger et al. 2004), a speckle- nulling experiment hosted at the Jet 
Propulsion Laboratory, specifically designed to test TPF-C related technology. 

After reviewing the process of speckle formation to establish our notations (§2), we 
derive two speckle nulling methods in the case of small aberrations (§3). The speckle nulling 
phase is preceded by the measurement of the electric field in the image plane (§3.4). The 
performance of both methods are then evaluated with one- and two-dimensional simulations 
(§4), first with white speckle noise (§4.1), then with non-white speckle noise (§4.2). Various 
effects and instrumental noises are considered in §5. Finally, we conclude and discuss some 
future work (§6). 

2. Speckle formation 

This paper is written in the framework of Fourier optics considering a single wavelength, 
knowing that a more sophisticated theory (scalar or vectorial) in polychromatic light will 
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eventually be needed. Fourier transforms (FTs) are signaled by a hat. 

Let us consider a simple telescope with an entrance pupil V . In the pupil plane, we use 
the reduced coordinates (u,v) = (x/X,y/X), where (x,y) are distances in meters and A is 
the wavelength. We define the pupil function by 

P(u,v)^{l [t ^ eV ' (1) 
v ; 10 otherwise. v ; 



Even in space, i.e. when not observing through a turbulent medium like the atmosphere, 
the optical train of the telescope is affected by phase and amplitude aberrations. Phase 
aberrations are wavefront corrugations that typically originate in mirror roughness caused by 
imperfect polishing, while amplitude aberrations are typically the result of a heterogeneous 
transmission or reflectivity. Moreover, Fresnel propagation turns phase aberrations into 
amplitude aberrations, and the reverse (e.g. Guyon 2005). Regardless of where they originate 
physically, all phase and amplitude aberrations can be represented by a complex aberration 
function in a re-imaged pupil plane, so that the aberrated pupil function is now Pe 1 ^. 

The electric field associated with an incident plane wave of amplitude unity is then 

E(u,v) = P(u,v)e i ^ u ' v \ (2) 



Exoplanet detection requires that we work in a regime where aberrations are reduced to 
a small fraction of the wavelength. Once in this regime, we can replace by its first order 
expansion 1 + i(f> (we will discuss in §5.2 the validity of this approximation). The electric 
field in the image plane being the FT of (2), we get 

E(a,/3)=P{a,l3) + iP<j>(a,P), (3) 

where (a, f3) are angular coordinates in the image plane. 

The physical picture is as follows. The first term (P) is the direct image of the star. 
The second term (P0) is the field of speckles surrounding the central star image, where each 
speckle is generated by the equivalent of first-order scattering from one of the sinusoidal 
components of the complex aberration 0. Each speckle is essentially a ghost of the central 
PSF. 

In the remainder of this paper, we focus on means to measure and correct the speckles 
in a coronagraphic image. Following Malbet, Yu, & Shao (1995) we will leave out the 
unaberrated PSF term by assuming that it is was canceled out by a coronagraph of some 
sort (see Quirrenbach (2005) for a review on coronagraphs) . Thus we clearly separate the 
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gain in contrast that can be obtained by reducing the diffracted light with the coronagraph 
on one hand, and by fighting the scattered light with the speckle nulling technique on the 
other hand. 



3. Speckle nulling theory 

The purpose of speckle nulling is to reduce the speckle noise in a central region of the 
image plane. This region, the dark hole, then becomes dark enough to enable the detection 
of companions much fainter than the original speckles. Speckle nulling is achieved by way 
of a servo system that has a deformable mirror as actuator. Because our sensing method 
requires DM actuation and is better understood with the knowledge of the command control 
theory, we first model the deformable mirror (§3.1), then present two algorithms for the 
command control (§3.2 & 3.3), and conclude with the sensing method (§3.4). 

3.1. Deformable mirror 

The deformable mirror (DM) in Trauger et al. (2003) consists of a continuous facesheet 
supported by iV x N actuators arranged in a square pattern of constant spacing. This 
DM format is well adapted to either square or circular pupils, the only pupil shapes that 
we consider in this paper 1 . We assume that the DM is physically located in a plane that is 
conjugate to the entrance pupil. However, what we call DM in the following is the projection 
of this real DM in the entrance pupil plane. The projected spacing between actuators is 
denoted by d. We assume that the optical magnification is such that the DM projected size 
is matched to the entrance pupil, i.e. Nd = D, where D is either the pupil side length 
or its diameter. The DM surface deformation in response to the actuation of actuator 
(k, I) G {0 . . . N— l} 2 is described by an influence function, denoted by The total phase 
change introduced by the DM (DM phase function) is 

N-l N-l 

ip(u,v) = ^2^2a k if M (u,v), (4) 

fc=0 1=0 

where au are actuator strokes (measured in radians). Note that contrary to the complex 
aberration function 0, the DM phase function is purely real. 



1 Two square DMs can be assembled to accommodate an elliptical pupil such as the one envisioned for 
TPF-C. 
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With an ideal coronagraph and a DM, the image-plane electric field formerly given by 
(3) becomes 

E'(a,P) = iP(p(a,P)+iPlp(a,P). (5) 

In the next two sections, we explore two approaches for speckle nulling. In §3.2, we 
begin naively by trying to cancel E' . Because there is a maximum spatial frequency that 
the DM can correct for, the DH has necessarily a limited extension. Any energy at higher 
spatial frequencies will be aliased in the DH and limit its depth. Therefore, the DM cannot 
be driven to cancel E', unless Pep is equal to zero outside the DH (i.e. unless there are 
already no speckles outside the DH). With this in mind, we start over in §3.3 with the idea 
that speckle nulling is better approached by minimizing the field energy. 



3.2. Speckle field nulling 

The speckle field nulling approach consists in trying to null out E' in the DH region 
(H), meaning we seek a solution to the equation 

V(a,/3)eH, P${a, 0) + ffria, P) = 0, (6) 

although, as we shall show, this equation has no exact solution unless P(p happens to be a 
band-limited function within the controllable band of the DM. 

By replacing ip with its expression (4), we obtain 

JV-1JV-1 

v(a, p)en, akl p /«( a > P) = ft- ( ? ) 

A:=0 1=0 

We recognize in (7) a linear system in the an that could be solved using various techniques 
such as singular value decomposition (SVD; Press et al. 2002, §2.6). Although general, this 
solution does not provide much insight in the problem of speckle nulling. For this reason, 
let us examine now a different solution, less general but with more explanatory power. We 
will comment on the use of SVD at the end of this section. 

We consider a square pupil. In this case, all DM actuators receive light and the pupil 
function has no limiting effect on the DM phase function, i.e. Pip = ip. Moreover, we assume 
all influence functions to be identical in shape, and write fki(u, v) = f(u — kf, v — Under 
these hypotheses, 

jv-ijv-i / d d\ 
Pip(u, v) = f(u, v)*J2J2 aklS [ u - k y v - l \) ' ( 8 ) 

k=0 1=0 ^ ' 



-7- 



where 5 is Dirac's bidimensional distribution, and * denotes the convolution. 
Substituting by the FT of (8) in (6) yields 

Ar ~ 1Ar - 1 ul( a\ 

V(a,0 e W, J2J2a kl e-^^ = -^l. (9) 
k=o i=o f( a 'P) 

We recognize in the left-hand side of (9) a truncated Fourier series. If we choose the a k i to 
be the first N 2 Fourier coefficients of —P(f>/f, i.e. 

«« = ^// 2 -S^f e^^aaap, (10) 

then according to Fourier theory, we minimize the mean-square error between both sides 
of the equation (see e.g. Hsu 1967, §1.5). This error cannot be reduced to zero unless the 
Fourier coefficients of —P(j)/f happen to vanish for k, I < and k,l > N. At this point, we 
have reached the important conclusion that perfect speckle cancellation cannot be achieved 
with a finite-size DM unless the wavefront aberrations are band- limited. Moreover, we can 
assert that the maximum DH extension is the square domain H = [— ^, ^] = hf ^1 ■ 

Solution (10) is physically acceptable only if the Fourier coefficients are real numbers, 
which means mathematically that Pep/ f should be Hermitian 2 . If there are phase aberrations 
only, P<p is real, P<f)/ f is Hermitian, and the a k i are real. This is no longer true if there 
are amplitude aberrations as well, reflecting the fact that the DM alone cannot correct both 
phase and amplitude aberrations in Ti. However, by considering the Hermitian function 
that is equal to P(f>/ f in one half of the DH, say H + = [0, ^] x [— ^, ^], we obtain the real 
coefficients ^ 

a kl =4d 2 [[ - P f^^ C o S [^(fca + dad/3, (11) 
JJh+ f(a,P) L A 

that correct both amplitude and phase aberrations in H + . As we have ^ = f A, the DH 
has a size of iV x N resolution elements (resels) with phase aberrations only, and of y x iV 
resels with phase and amplitude aberrations. Therefore, a DM can correct both amplitude 
and phase aberrations in the image plane, albeit in a region that is either the left, right, 
upper, or lower half of the phase- corrected region. 

As Malbet, Yu, & Shao (1995) pointed out, let us remind the reader that A i s equal to 
the Nyquist frequency for a sampling interval 4. Therefore, we find that maximum extension 



2 A function / is said to be Hermitian if V(x,y), f(x,y) = f*(—x,—y). The FT of a real function is 
Hermitian and vice versa. 
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for the DH corresponds to the range where the sampling theorem applies to the wavefront 
at the DM actuator scale. Indeed, taking the inverse FT of (9) leads to the wavefront 
reconstruction formula 

n-in-i / H H\ 
P<f>(u, v) = - £ «« / ( « - kj,v - lj J . (12) 

k=0 1=0 ^ ' 

Again, this reconstruction cannot be perfect unless the spectrum of P0 is contained in Ti. 
Note that because / is generally not a flat function (as it would be the case if influence 
functions were for instance 2D sine functions), actuator strokes are not equal to the negative 
of wavefront values sampled at the actuator locations. 

Our Fourier solution was derived by assuming that (a) all influence functions are iden- 
tical in shape, and (b) that the pupil has a square shape. Hypothesis (a) appears to be 
reasonable at least for the DM in use on the HCIT (Joseph Green, personal communica- 
tion), but this remains to be precisely measured. If hypothesis (b) is relaxed then (i) some 
actuators do not receive any light and play no role, so there are effectively fewer terms in 
the summation in (9), and (ii) the fact that influence functions on the pupil boundary are 
only partly illuminated is ignored. 

Now that we have two methods to solve (6), Fourier expansion and SVD, let us com- 
pare their solutions. We deal here with functions belonging to the Hilbert space of square 
integrable functions / : Ti — > C. This space has < /, g >= ff n f g* for dot product, and 

ll/H = yj JJ n |/| 2 for norm. As mentioned earlier, Fourier expansion minimizes the mean- 
square error between both sides of (9), i.e. ||(P0 + P0)//|| 2 . By contrast, SVD has the 
built-in property of minimizing the norm of the residuals of (7), i.e. \\P<f> + Pip\\- In other 
words, SVD minimizes ||-E'|| 2 , the speckle field energy, which seems more satisfactory from 
a physical point of view. To find out what is best, we have performed one-dimensional nu- 
merical simulations. It turns out that SVD yields dark holes 50 % deeper (median value) 
than Fourier expansion. In addition, SVD does not require all influence functions to have 
the same shape. 

However, considering four detector pixels per resel in two dimensions (critical sampling), 
SVD would require us to manipulate matrices as large as A^ 2 x 4A^ 2 (or even A^ 2 x 8 N 2 when 
real and imaginary parts are separated). Such matrices would occupy 537 MB of memory 
space for 64x64 actuators and single-precision floating-point numbers. By contrast, Fourier 
expansion would be straightforwardly obtained with FFTs of 2A" x 2A" arrays at critical 
sampling, but again at the cost of a 50 % shallower dark hole and a strong hypothesis on the 
influence functions. 

In the next section, we seek to find a computationally less intensive solution that still 
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minimizes the speckle energy in the dark hole, but does not require any hypothesis on the 
influence functions. 



3.3. Speckle energy minimization 

Let us start with the idea that the best solution is defined as the one minimizing the 
total energy of the speckle field in the DH. For the sake of simplicity, we assume once again 
a square pupil, but not necessarily a common shape for the influence functions. The total 
energy in the speckle field reads 

£ = J J +$(a,[3)\ 2 dad[3 = < P0 + P0 + ? >, (13) 

using the same notation as in §3.2. 

Given that dip/daki = fki, the energy is minimized when 

v(M)e{o...iV-i} 2 , ^ = o k(<P0 + ?,/«>)=o, (14) 

where 3ft stands for the real part. Note that this is less demanding than (6), as (6) implies 
(14) but the reverse is not true. 

Using the definition (4) for ip and realizing that < f nm , fki > is a real number 3 , we get 
finally 

N-l N-l 

V(fc, I) e {0 . . . N-l} 2 , E a ™ < f™> /« > = (< P<j>, /« >) . (15) 

n=0 m=0 

As in (7), we find a system that is linear in the actuator strokes. By replacing double 
indices with single ones, e.g. (k, I) becomes s — k N + 1, (15) can be solved in matrix format 
by inverting a iV 2 x iV 2 real matrix. This is already an improvement with respect to the 
iV 2 x 4iV 2 complex matrix required by SVD in the previous section. 

It appears that the same solution can be obtained with a much less demanding NxN 
matrix inversion, provided two-dimensional influence functions can be written as the ten- 
sor product of two one-dimensional functions (separation of variables), i.e. fki{u,v) = 
gk(u)gi(v). This would be the case for box functions or bidimensional Gaussians, and is 



3 This property stems from the Hermitian character of fki together with the symmetry of Tl. 
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good at the 5 % level for the DM in use on the HCIT. This property also holds in the image 
plane since the FT of the previous equation yields fki(a,/3) = gk(&) gi(/3). 

By separating variables, (15) becomes 

N-l N-l 

V{k,l) e {0...N-1} 2 , J2<9n,9k> J2 a ™<9rn, 9i > = P<f>, hi >) ■ (16) 

n=0 m=0 

As the left-hand side happens to be the product of three NxN matrices, (16) can be rewritten 
as an equality between square matrices. 

{G k i = < g k ,gi > 
A M = a kl (17) 
= (< P<t>J k i >). 

For square-box and actual HCIT influence functions, numerical calculations show that G is 
diagonally dominant 4 and therefore invertible by regular Gaussian elimination. The solution 
to (17) is then 

A = G- 1 <$>G- 1 . (18) 

Note that G^ 1 can be precomputed and stored, so that computing the strokes effectively 
requires only two matrix multiplications. As shown in appendix A, an equivalent result can 
be obtained by working with pupil plane quantities. 

As for the field nulling approach, correcting amplitude errors as well implies restricting 
the dark hole to either H + = [0, f A] x [-f A |A] or H - = [_A A ] x [-f A AAj. To 
account for amplitude errors and keep the formalism we have presented so far, it is sufficient 

to replace P<p by a function equal to P<f>(a,/3) in either 7i + or 7i~ (depending on the half 

- — - * 

where one wishes to create the dark hole), and equal to P<p (—a,— p) in the other half 
(Hermitian symmetry). Because its FT is Hermitian, the new aberration function in the 
pupil plane is real, and thus the algorithm processes amplitude and phase errors at the same 
time as if there were phase errors only. 

Let us derive the residual total energy in the DH after the correction has been applied. 
Starting from definition (13) and rewriting condition (14) as 9R(< P<j) + ip, ip >) — 0, we find 

£min = <P0,P0>-<V^>. (19) 

The former term is the initial speckle energy in the DH, while the latter is the speckle energy 
decrease gained with the DM. Mathematically, V^-min measures the distance (according to 



A matrix A = [oy] is said to be diagonally dominant if Mi, \au\ > J2j^i \ a ij 
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the norm we have defined) between the speckle field and its approximation with the DM 
inside the DH. Because there is no exact solution to (6), the residual energy cannot be made 
equal to zero in 7i + or . However, the energy approach offers an additional degree of 
freedom: by reducing concentrically the domain over which the energy is minimized, the 
speckle energy can be further decreased (see §4). 



3.4. Speckle field measurement 

So far, our speckle nulling theory has presupposed the knowledge of the speckle field 
P(f>, or equivalently of the phase and amplitude aberrations across the pupil, embodied in the 
complex phase function P<fr. In this section, we show how the speckle field can be measured 
directly in the image plane. As the detector measures an intensity, a single image yields only 
the modulus of the speckle field. The phase of the speckle field can be retrieved by perturbing 
the phase function P(f> in a controlled way, and by recording the corresponding images, a 
process analogous to phase diversity (e.g. Lofdahl & Scharmer 1994). In our system, the 
DM provides the natural means for creating this controlled perturbation. 

As we will see, exactly three images obtained with well chosen DM settings provide 
enough information to measure P(p. Let us call image the original image recorded with a 
setting tp , whereas images 1 and 2 are recorded with settings ip + 5ipi and ip + 5ifj 2 - 

To be general, we consider in the field of view the presence of an exoplanet and an 
exozodiacal cloud (exozodi for short), in addition to the star itself. The electric fields of 
these objects are incoherent with that of the star, so their intensities should be added to 
the star's intensity. Because they are much fainter than the star, the speckles they produce 
are negligible with respect to the star speckles, and their intensities can be considered as 
independent of and ip. The total intensity of every image pixel (a, (3) takes then the 
successive values 

J = |P0 + ^ O | 2 + /p + /z 

/ 1 = |P0 + ^o + §S| 2 + / P + /z (20) 

h = \P<t> + $ + fohl 2 + / p + / z , 

where I p and I z are the exoplanet and exozodi intensities, respectively. 
System (20) can be reduced to the linear system 

(go* (P0 + ft) + gj W + ft)* =h-h- ig^i 2 (21) 

(<^ 2 )* {Pct> + ft) + 5ft w + ft)* =h-h- i^Si 2 , 

where the exponent * denotes the complex conjugate. Notice how the exoplanet and exozodi 
intensities have disappeared from the equations, demonstrating that faint objects do not 
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affect the measurement process of stellar speckles. However, note that because of quantum 
noise, the planet detection can still be problematic if the exozodi is much brighter than the 
planet. 



Now, system (21) admits a unique solution if its determinant, 

A = (^)*^-^(^)*, 

is not zero, that is to say if 

|#0i| \Sip2\ sin arg(#0 2 ) - arg(#0i) ^ 0. 



(22) 



(23) 



Condition (23) tells us that the DM setting changes, Sipi and 5ip2, should modify the 
speckles differently in any given pixel, otherwise not enough information is secured to measure 
unambiguously P<p in this pixel. It should be expected for this method to work practically 
that the magnitude of the speckle modification be greater than the photon noise level. 

We have not yet found a rigorous derivation of the optimum values for the amplitude 
| | and 1 8ip2 1) but a heuristic argument suggests to us that the optimum perturbations 
may be that I\ ~ Iq and I 2 ~ Jo- That is to say, the DM-induced speckle instensity pattern, 
taken by itself, should be approximately the same as the original speckle intensity pattern. 
Thus at each pixel we choose \5ipi\ ~ |#0 2 | ~ V^o, with the caveat that neither should be 
zero to keep (23) valid. 

The phase of 5ipi does not matter, but the phase difference between Sipi and 5ip 2 should 
be made as close to f as possible to keep A from zero. Practically, this can be realized as 
follows: 

1. Compute 5tpi stroke changes from (11) or (18) by replacing P<p by \H^e t6 , where 9 is 
a random phase; 

2. Compute 5ip 2 stroke changes from (11) or (18) by replacing P<j) by e ! i 

Now that we have made sure that A ^ 0, we derive finally 

Si (h -Iq- I W) -jMh-Jo- I W) T f0A . 
P(p = — O . (24) 

Equation (24) shows that the initially unknown speckle field (P0) can be experimentally 
measured in just three exposures taken under identical circumstances but with different 
shapes imposed on the DM. 
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4. Speckle nulling simulations 



4.1. White speckle noise 



In this section, we perform one- and two-dimensional simulations for the theoretical 
case of white speckle noise caused by phase aberrations only. The DM has 64 actuators and 
top-hat influence functions. Smoother influence functions have been tested and do not lead 
to qualitatively different results. A simulation with actual HCIT influence functions will be 
presented in the next section. The simulated portion of the pupil plane is made twice as big 
as the pupil by zero padding, so that every element of resolution in the image plane would 
be sampled by two detector pixels. This corresponds to the realistic case of a photon-starved 
exoplanet detection where read-out noise must be minimized. 



Figure 1 shows a complete one-dimensional simulation including speckle field measure- 
ment (§3.4) and speckle nulling with field nulling (§3.2) and energy minimization (§3.3). 
The standard deviation of the phase aberrations is set to A/1000. Intensities are scaled with 
respect to the maximum of the star PSF in the absence of a coronagraph. Ideal conditions 
are assumed: no photon noise, noiseless detector, and perfect precision in the control of DM 
actuators. Under these conditions, the speckle field is perfectly estimated, and the mean in- 
tensity in the DH is 5.8 x 10~ n with field nulling and 1.4 x 10~ n with energy minimization, 
i.e. about 1500 and 6500 times lower than the mean intensity outside the DH, respectively. 
Repeated simulations with different noise sequences show that energy minimization performs 
always better than field nulling by a factor of a few. Field nulling solved with SVD yields 
the same numerical solution as energy minimization (they differ by the last digit only), in 
agreement with the idea that they both minimize speckle energy. 



In the one-dimensional case, it is easy to predict roughly the shape and the depth of 
the DH. The function P0 + ip is band-limited since the pupil has a finite size. As the pupil 
linear dimension is D/X, the maximum spatial frequency of Pip + ip is D/2X. Let us apply 
the sampling theorem at the Nyquist sampling frequency D/X, and write 



4-1.1. One- dimensional simulations 



4-1.2. Dark hole depth estimate in one dimension 




n=—oo 



(25) 
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where the subscript n denotes the function value for a = n^. 
Substituting a by n-^ and d by leads to 

N-l 

P0 n + ? n = P0 n +/ ri ^ a fe e-^ . (26) 



fc=0 

The field nulling equation (6) takes here the discrete form 

i.e. actuator strokes are computed thanks to an inverse FFT. 
Let us now turn to the residual speckle field 



{P<i> + $)(a)= ^>nsinc(^--nj + ^ P0 n 

n=— oo ^ ' n=N 



sinc(^-n) (28) 



Because the sine function decreases rapidly with a, the terms flanking the DH (n — — 1 and 
n = N) should by themselves give the order of magnitude of the residual speckle field in 
the DH. In case of phase aberrations only and white noise, we have |P0_ 1 | 2 \P(f) N \ 2 ~ Iq, 
where Jo is the mean intensity in the image plane prior to the DH creation. Therefore, a 
crude estimate of the intensity profile in the DH should be 

2 



'DH 



a) 



. aD \ . aD Ar 
sine ( — - — h 1 I + sine I — JM 



(29) 



We have superimposed this approximation as a thick line in Fig. 1. In this case the match 
is remarkable, but more simulations show that it is generally good within a factor of 10 only. 
Nevertheless, it demonstrates that the DH depth depends critically on the residual speckle 
field at its edges, hence on the decreasing rate of the complex aberration spectrum with 
spatial frequency. In that respect, a white spectrum is certainly the worst case. Equation 
(29) further indicates that the DH depth depends on the number of actuators: as iV is 
increased, the DH widens and gets deeper. With 8, 16, 32, and 64 actuators, (29) predicts 
To/Tdh" to reach about 100, 300, 1000, and 4500. 



4-1.3. Dark hole depth vs. search area 



As mentioned in §3.3, speckle nulling by energy minimization can be performed in 
a region narrower than the maximum DH. Figure 2 illustrates this point: by reducing the 
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search area from 64 to 44 resels (31 % reduction), the DH floor was decreased from 1.4 x 1CT 11 
to 2.7 x 10~ 15 , i.e. a gain of about 5200 in contrast (further reducing the search area 
does not bring any significant gain). By giving up search space, one frees the degrees of 
freedom corresponding to the highest spatial frequency components on the DM pattern. 
These can be used to improve the DH depth at lower spatial frequency because of the PSF 
angular extension (this is essentially the same reason why high spatial frequency speckles 
limit the DH depth). As the search space is reduced, the leverage of these highest spatial 
frequency components decreases (PSF wings falling off). The energy minimization algorithm 
compensates by putting more energy at high frequency (see lower panel in Fig. 2), which 
produces increasingly oscillatory DM patterns (see top panel in Fig. 2) and increasingly 
brighter spots in the image (around ±32-^ and ±96-^ in bottom panel of Fig. 2). Thus 
the trade-off range might be limited in practice by the maximum actuator stroke (currently 
0.6 /iin on the HCIT), and/or by the detector's dynamic range. 

In two-dimensions, the trade-off limits are well illustrated by the following example: 
considering a 64x64 DM and a random wavefront, we find that the DH floor can be decreased 
from 2.4 x 10~ 12 to 1.4 x 10~ 13 (a factor of 17) if the search area is reduced from 64x64 to 
60x60 resels (12% reduction in area). This implies a maximum actuator stroke of 10 nm 
and a detector dynamic range of 10 6 . A further reduction to 58x58 resels does not feature 
a lower DH floor (2.1 x 10~ 13 ), and would imply a maximum actuator stroke of 10 fxm and 
a detector dynamic range of 10 10 . In this case, the leverage of the additionally freed high- 
spatial frequency components is so weak that the algorithm starts diverging. 

Two-dimensional simulations with phase and amplitude aberrations 

In Figs. 3-4, we show an example of two-dimensional speckle nulling with phase and 
amplitude aberrations for a square pupil. To reflect the fact that phase aberrations dominate 
amplitude aberrations in real experiments (see Trauger et al. 2004), the rms amplitude of 
amplitude aberrations is made ten times smaller than that of phase aberrations (the choice of 
a factor ten is arbitrary). The DH is split into two regions: in the right one (H + ), amplitude 
and phase aberrations are corrected, whereas in the left one (H~), phase aberrations are 
corrected and amplitude aberrations are made worse by a factor of four in intensity. 
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4.2. Realistic speckle noise 

4-2.1. Power spectral density of phase aberrations 

With the 3.5x8-m TPF-C primary mirror in mind, we have studied the phase aberration 
map of an actual 8-m mirror: the primary mirror of Antu, the first 8.2-m unit telescope of 
ESO's Very Large Telescope (VLT). This phase map 5 was obtained with the active optics 
system on, and is characteristic of zonal errors (aberrations which cannot be fitted by low- 
order Zernike-type polynomials). It can be seen in Fig. 5 that the azimuthally averaged 
power spectral density (PSD) of such a map is well represented by 

PSD 

PSD <" = T+Jpjpcf' (30) 

where p = \J a 2 + (3 2 . Values for PSDo, p c and x are listed in Table 1. For comparison, 
the same treatment has been applied to the Hubble Space Telescope (HST) zonal error map 
from Krist & Burrows (1995). 

We conclude from this study that a realistic phase aberration PSD for an 8-m mirror 
decreases as the third power of the spatial frequency. The standard deviation of the VLT 
phase map is 20.9 nm (18.5 nm for HST). The square root of the power of phase aberrations 
in the 0.5-4 m -1 spatial frequency range (4-32 X/D for an 8-m mirror) is 19.4 nm, i.e. about 
A/25 at 500 nm, clearly not in the validity domain of our linear approximation. 



4-2.2. One-dimensional simulation 

Figure 6 shows a simulation performed in the same conditions as Fig. 1, but with a 
VLT-like PSD. The PSD is scaled so that the standard deviation of phase aberrations is 
equal to A/1000. The average DH floor is now 5.3 x 10~ 12 , six orders of magnitude below the 
intensity peak in the original image! In agreement with §4.1, we find that the DH's depth 
depends critically on the magnitude of the speckle field at the edge of the DH, hence on the 
decrease of the phase aberration PSD with spatial frequency. 



5 It can be found by courtesy of ESO at http://www.eso.org/projects/vlt/unit-tel/mlunit.html. 
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4-2.3. Two-dimensional simulation 

For the two-dimensional simulation in Figs. 7-8, we have kept the original VLT phase 
map and circular pupil, but scaled the standard deviation of phase aberrations to A/ 1000. 
In addition, we have used the actual HCIT influence functions from Trauger et al. (2003). 
The average DH floor is then 5.9 x 10~ 12 with field nulling (case shown), and 7.1 x 10~ n 
with energy minimization. The worse performance of the second method reflects the cost of 
the variable separation hypothesis, only accurate to within 5 % for the HCIT. Note that the 
DH retains its square shape with a circular pupil, as the DH shape is fixed by the actuator 
grid geometry on the DM (a square grid of constant spacing in our case). 



5. Discussion 

5.1. Quantum and read-out noise 

In §4, we presented noise-free simulations. To give an idea of the effect of quantum and 
read-out noises, let us consider a sun-like star at 10 pc observed by a 3.5x8 m space telescope 
with a 5 % overall efficiency In a 100 nm bandwith centered at 600 nm, the telescope collects 
about 2 x 10 12 photo-electrons in one-hour exposures. Considering the quantum noise, ale" 
read-noise and ignoring chromatic effects, simulations of sequences of four one-hour exposures 
show that the average DH floor in Fig. 1 would jump from 1.4 x 10 -11 to 2.7 x 10~ 10 , whereas 
the average DH floor in Fig. 6 would jump from 5.2 x 10~ 12 to 3.2 x 10~ n . 



5.2. Validity of the linear approximation 

In practice, our speckle nulling process will work as stated provided Eq. (3) holds, 
that is to say if \P(f> + ip\ ^> ||P0 2 |. If c is the improvement in contrast with respect to the 
speckle floor and the standard deviation of wavefront aberrations in radians, this condition 
translates into a^/y/c^ or <C \/2/c. In terms of optical path difference, the 

standard deviation should then be much less than \/(ix\f2c) = A/140 for c = 10 3 . This is 
why we considered A/1000 rms wavefronts in our simulations. As the wavefront will probably 
not be of this quality at the start, the speckle nulling method presented here is intended to 
be used in the course of observations, after a first phase where the bulk of the aberrations 
have been taken out. 

When the linear approximation breaks down, three images with different DM settings 
still provide enough information about the aberrations, so that a DH could be created thanks 
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to a global non-linear analysis of these images (Borde, Traub, & Trauger 2004). Malbet, Yu, 
& Shao (1995) also explored non-linear solutions, but with many more iterations 20). 

5.3. Real coronagraphs 

Dwelling on the validity of Eq. (3), real coronagraphs would not only remove the direct 
image of the star (P), they would also modify the speckle field (P0) and the DM phase 
function (Pip). This can be easily incorporated in the theory. A more delicate point is 
that real coronagraphs are not translation-invariant systems. As a consequence, effective 
influence functions as seen from behind the coronagraph will vary over the pupil. For image- 
plane coronagraphs with band-limited sine masks (Kuchner & Traub 2002, §4), we estimate 
this variation to be of the order of 10%, assuming e = 0.1 and 64 actuators. Only energy 
minimization, not field nulling (unless solved with SVD), can accomodate this effect. 

5.4. Actuator stroke precision 

What about the precision at which actuators should be controlled? As a consequence of 
the linearity of (15), the DH depth depends quadratically on the precision on the actuator 
strokes. We deduce - and this is confirmed by simulations - that a four orders of magnitude 
deep DH can only be obtained if the strokes are controlled at a 1 % precision, i.e. 6 pm rms 
with A/ 1000 aberrations at 600 nm. This precision corresponds to the current resolution of 
the actuator drivers on the HCIT. 



5.5. Instrumental stability 

Regarding instrumental stability, we assumed that the instrument would remain per- 
fectly stable during the four-step process. However, despite the foreseen thermal and me- 
chanical controls of the spacecraft, very slow drifts during the few hours of single exposures 
should be expected. Therefore we intend to study in a subsequent paper how to incorporate 
a model of the drifts in our method. The exact parameters of this model would be derived 
from a learning phase after the launch of the spacecraft. 
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5.6. Chromaticity 

We have not considered the effect of chromaticity. Let us point out that phase aberra- 
tions due to mirror surface errors scale with wavelength, so the correction derived from one 
wavelength would apply to all wavelengths. This is unfortunately not the case for ampli- 
tude aberrations. Although these are weaker than phase aberrations, a degradation of the 
correction should be expected in polychromatic light. Moreover, polychromatic wavefront 
sensing will require a revised theory as speckles will move out radially in proportion to the 
wavelength. 

6. Conclusion and future work 

In this paper, we presented two techniques to optimally null out speckles in the central 
field of an image behind an ideal coronagraph in space. The measurement phase necessi- 
tates only three images, the fourth image being fully corrected. Depending on the number 
of actuators and the desired search area, the gain in contrast can reach several orders of 
magnitude. 

These techniques are intended to work in a low aberration regime, such as in the course 
of observations after an initial correction phase. They are primarily meant to be used in 
space but could be implemented in a second-stage AO system on ground-based telescopes. 
Out of these two methods, the speckle energy minimization approach seems to be the more 
powerful and flexible: (i) it offers the possibility to trade off some search area against an 
improved contrast, and (ii) it can accomodate influence function variations over the pupil 
(necessary with real coronagraphs) . If influence functions feature the required symmetry 
(variable separation), it is computationally very effective, but is otherwise still better than 
SVD. 

Since the principles underlying these speckle nulling techniques are general, it should 
be possible to use them in conjunction with most coronagraph designs, including those with 
band-limited masks (Kuchner & Traub 2002), pupil-mapping (Guyon et al. 2005; Vanderbei 
& Traub 2005), and shaped pupils (Kasdin et al. 2003). It is our intent to complete our work 
by integrating in our simulations models of these coronagraphs, and to carry out experiments 
with the HCIT. 

In addition, we will seek to incorporate in the measurement theory a linear model for 
the evolution of aberrations, and we will work toward a theory accommodating the spectral 
bandwidth needed for the detection and spectroscopy of terrestrial planets. 
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A. Formulation of energy minimization in the pupil plane 

In this appendix, we show that energy minimization can be formulated in the pupil 
plane as well. Note that no measurement takes place in the pupil plane: the aberration 
function P<p is obtained by the inverse FT of P<p which is still measured as described in 
§3.4. Although we present here a solution for phase aberrations, amplitude aberrations can 
be corrected in half of the domain without changing the formalism, exactly as explained in 
§3.3. 

By replacing E with its expression as a FT, the dark hole energy reads 
£ = jj E(u,v)e- t27T{ua+v(S) dudv J J E*(u',v')e t27T{u ' a+v,fS) du'dv^dad[3. (Al) 

Now we invert the integration order and integrate over Ti to get 

£ = -LJJ E(u,v) J J E*{u',v')h{u' -u)h{v' - u) dw du d«W, (A2) 
where h(x — y) = sinc[A(a; — y)/2d\. The energy is minimized when 

V(Jfe,Z), // fki(u,v) [[ h(u',v')+il>(u',v') \ h(u' -u)h(v' -v)dudvdu'dv' = 0. (A3) 



In the next steps, we first replace ip with (4), then fki{u,v) with its tensor product 
9k( u ) 9i( v ) m order to integrate separately the variables (u,u') and (v,v r ). The final result 
reads 

N-1 N-1 

V(fc, I) e {0 . . . N-l} 2 , G >™ Yl a ™ Gml = ® kl ( A4 ) 

n=0 m=0 

where = J J g t (x) g^y) h(x - y) dxdy 

and $m = J J gk{u) gi(v) J J 4>(u', v') h{u — u) h(v' — v) du dv du'dv'. 
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System (A4) has a form identical to (17) and can be solved with the same technique. 
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Fig. 1. — Full one-dimensional speckle nulling simulation for a one-dimensional pupil with 
64 actuators. Panels on the left show phase aberrations as dotted lines and their low spatial 
frequency approximations at the actuator scale as solid lines (negative of the patterns applied 
to the DM). Panels on the right show the corresponding images. The full algorithm is a four- 
step process: (a) original speckles are measured with the current DM shape (taken here to 
be flat); (b) and (c) the speckles are modified by driving the DM. At this point enough 
information has been gathered to deduce phase and amplitude aberrations of the wavefront; 
(d) low spatial frequency aberrations are corrected with the DM, canceling out speckles in 
the center of the image. Dark holes produced by the field nulling and energy minimization 
approaches (§3.2 & 3.3) do not differ noticeably on the figure (thin solid line). A rough 
estimate based on the mean intensity prior to the dark hole creation (§4.1) is superimposed 
as a thick solid line. 
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Fig. 2. — The energy minimization algorithm makes it possible to push down the dark hole 
floor at the cost of some search area: the original dark hole appears in thick solid line, while 
the deeper and narrower dark hole appears in thin solid line. In this particular example 
with phase aberrations only, the average floor is decreased from 1.4 x 10~ n to 2.7 x 1CT 15 (a 
factor of about 5200 in contrast) by reducing the dark hole size from QAX/D to 44A/-D (31 % 
reduction). As explained in the text, this trade-off is obtained by increasing the amplitude 
of the highest spatial frequencies on the DM, making the dark hole's rim brighter at the 
same time. 



-25 - 



Dark hole (64 x64 actuators) 




ailID) 

Fig. 3. — Speckle nulling with the energy minimization algorithm (§3.3) for a square pupil 
in two dimensions. Grey levels code the logarithm of the intensity. Speckles are cleared 
from the central part of the image, creating a dark hole (or search area) suitable for the 
detection of faint companions. Phase aberrations are corrected in the full dark hole, while 
amplitude aberrations are corrected in the right part only and made worse by a factor of 
four in intensity in the left part. Thus the difference in intensity between the two sides of 
the dark hole gives a measure of the wavefront amplitude errors. In this simulation, the 
standard deviations of phase and amplitude aberrations are A/10 3 and A/10 4 , respectively. 
The dark hole shape is a result of the actuator grid geometry on the DM, not of the pupil 
shape. 
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Average cut through the dark hole 
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Fig. 4. — Speckle nulling with the energy minimization algorithm (§3.3) in two dimensions. 
The solid curve is an average of Fig. 3 intensity over (3. The dotted line represents the 
state prior to correction with the DM. Notice that the dark hole has a rim brighter than the 
background. 
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Azimuthally averaged PSD 
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Fig. 5. — Azimuthally averaged power spectral density (PSD) of phase aberrations for (i) 
the 8.2-m primary mirror of a VLT unit telescope (Antu), and (ii) the combination of the 
primary and secondary mirrors of Hubble. Both data sets appear in solid line, and are 
fitted with model (30) drawn in dotted line. Model parameters are listed in Table 1. The 
dashed lines indicate the boundaries for spatial frequencies leading to speckles in the 4-32 
X/D region of the image plane for an 8-m mirror. The lower- left and upper-right insets 
are the VLT phase map (courtesy of ESO) and HST phase map (Krist & Burrows 1995), 
respectively. 
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Pupil plane (64 actuators) 
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Fig. 6. — One-dimensional speckle nulling as computed by energy minimization (§3.3) with a 
phase aberration PSD given by (30). Values for x and p c are the same as for the VLT; PSDo 
is such that the standard deviation of phase aberrations is the same as in Fig. 1 (A/1000). 
These realistic phase aberrations lead to deeper dark holes than in the hypothetical white- 
noise case of Fig. Id. 
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Fig. 7. — Speckle nulling with the field nulling algorithm (§3.2) with a VLT-like PSD instead 
of white noise. The standard deviation of phase aberrations is fixed to A/1000. Actual HCIT 
influence functions are used. The average DH floor is 5.9 x 10~ 12 . 
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Fig. 8. — Speckle nulling with the field nulling algorithm (§3.2) with a VLT-like PSD instead 
of white noise. The solid curve is an average of Fig. 7 intensity over (5. The dotted line 
represents the state prior to correction with the DM. 
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Table 1: Parameter values for the azimuthally averaged PSD model of VLT and HST phase 
maps. 





PSD (nm 2 m 2 ) 


Pc (m- 1 ) 


X 


HST 


2.2 


4.3 


2.9 


VLT 


720 


0.35 


3.1 



